Modelling of nematic liquid crystal display devices 
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A lattice Boltzmann scheme is presented which recovers the dynamics of nematic and chiral Hquid 
crystals; the method essentially gives solutions to the Qian-Sheng Q equations for the evolution 
of the velocity and tensor order-parameter fields. The resulting algorithm is able to include five 
independent Leslie viscosities, a Landau-deGennes free energy which introduces three or more elastic 
constants, a temperature dependent order parameter, surface anchoring and viscosity coefficients, 
flexo-electric and order electricity and chirality. When combined with a solver for the Maxwell 
equations associated with the electric field, the algorithm is able to provide a full 'device solver' for 
a liquid crystal display. Coupled lattice Boltzmann schemes are used to capture the evolution of 
the fast momentum and slow director motions in a computationally efficient way. The method is 
shown to give results in close agreement with analytical results for a number of validating examples. 
The use of the method is illustrated through the simulation of the motion of defects in a zenithal 
bistable liquid crystal device. 

PACS numbers: 61.30.-v,42.79.Kr 



I. INTRODUCTION 

Many of the next generation of liquid crystal (LC) dis- 
play devices use structured or patterned surfaces as an 
essential element of their design and function H, H, Q . 
The correct operation of these devices depends upon the 
formation and annihilation of defects in the orientational 
field of the nematic; the defects are usually intimately 
coupled to a surface. There is also current interest in 
the behaviour of LC's with embedded colloidal particles 
{eg P, 6]); the behaviour of these materials is frequently 
dependent upon the interaction of the defects associated 
with the colloidal particles. 

Experimentally it is difficult to obtain information 
about the spatial and temporal behaviour of the nematic 
order. Optical methods such as those of 0, ll| can, for 
example, give information about the director profiles on 
a relatively coarse time and space scale. However, in or- 
der to fully understand such systems it is necessary to 
be able to model the statics and dynamics of a nematic 
in the presence of complex boundaries and defects. The 
problem is compounded by the numerous materials pa- 
rameters needed to fully describe the properties of the LC 
and its interaction with any bounding surfaces; predictive 
modelling often requires a fairly complete description of 
the materials and this therefore necessitates the use of 
numerical methods to solve the associated equations. In 
this paper we present the details of one such numerical 
method and illustrate its use with a number of examples. 

LC's are complex fluids formed from anisometric 
molecules. These fluids can exhibit a range of mesophases 
with varying degrees of orientational and positional or- 
der of the molecules; in the nematic phase there is long 
range orientational order but no positional order. The 
orientational ordering is described at mesoscopic length 
scales by an order tensor, the Q-tensor (see Section . 
whose principal eigenvalue is related to the order param- 
eter and whose principal eigenvector deflnes the macro- 



scopic director fleld eg 0, Hol [Tlj . In many systems it is 
possible to assume that the order parameter is constant 
and the dynamics of the momentum and director are then 
described by the well established Ericksen-Leslie-Parodi 
(ELP) equations (e^ HJ). 

However, near to bounding walls and close to defects, 
the assumption of constant order parameter breaks down 
and the material may also exhibit biaxiality. There are 
signiflcant spatial gradients in the order tensor in such 
regions and the gradients have observable macroscopic 
consequences. For example, Q-tensor gradients lead to 
the flexo- and order-electric polarisation which is used to 
control the switching behaviour of some display devices 
by an applied electric field. Similarly, the dynamics of 
defects can only be correctly described within a theoret- 
ical framework which allows for variation in the order 
parameter. 

In such systems it is necessary to go beyond the ELP 
theory and adopt a model which describes the dynam- 
ics of the full Q-tensor. There are a number of deriva- 
tions of nematodynamics with variable order parameter 
{eg „1, 12, 13, 14]). Work by Sonnet et al provides the 
basis upon which the variety of schemes with a variable 
order parameter may be compared; it should be noted 
that in the limit that the order parameter becomes inde- 
pendent of time and position all the schemes must recover 
the ELP theory. In this work we adopt the Qian-Sheng 
formalism because it is straightforward to obtain the 
required material parameters from those of the equivalent 
ELP description. Exact solutions to either the ELP or 
Q-tensor theories are limited to a relatively small number 
of simplified cases and numerical methods are necessary 
for the more complex systems which form the focus of 
this work. 

A number of different approaches have been taken to 
finding numerical solutions to the equations for variable 
order parameter nemato-dynamics. Svensek p!^ and 
Fukuda ^3 use conventional methods to solve the as- 
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sociated partial differential equations. However, a num- 
ber of workers have adapted the lattice Boltzmann (LB) 
method (e^ [UlliiSlllllSl)- Care et al |3 developed 
the first methods in which LB was used to solve the ELP 
equations in a 2-D plane, and later [21j| enhanced the 
method to yield a two dimensional solution to the steady 
state Qian-Shcng equations; concurrently, Denniston et 
al {eg [20, 22, 23]) developed LB tensor methods for 
nematic liquid crystals based on the Beris-Edwards [l^ 
scheme. 

The LB method may be regarded simply as an alterna- 
tive method of solving a target set of macroscopic differ- 
ential equations. However, it is advantageous to regard 
it as a mesoscale method which allows additional physics 
to be included within the modelling; this is illustrated 
by the extension of the method to model the interface 
between an isotropic and nematic fluid 2lJ, a problem 
of direct relevance to modelling liquid crystal colloids. 
LB has the additional stability advantages of being able 
to incorporate complex boundary conditions more easily 
than conventional solvers and being straightforward to 
parallelise. 

In this paper we present an LB scheme which recovers 
the Qian-Sheng equations for nemato-dynamics. The ap- 
proach modifies the scheme presented in 21] by utilising 
a simple LBGK scheme {eg [2J]) for the collision term and 
introducing all the anisotropic behaviour through forcing 
terms. This moves away from the goal which was implicit 
in [IJ] of remaining as close as possible to the physical 
basis of nematodynamics by using an anisotropic collision 
operator. However, the overhead in numerical complexity 
made the scheme difficult to generalise to three dimen- 
sions, a restriction which does not apply to the method 
presented in this paper. 

The resulting algorithm includes five independent 
Leslie viscosities, a Landau-deGennes free energy which 
introduces three or more elastic constants, a tempera- 
ture dependent order parameter, surface anchoring and 
viscosity coefficients, flexo-electric and order electricity 
and chirality. The precise properties of the system, such 
as the number of elastic constants, is modified by the in- 
clusion or exclusion of terms in the free energy. When 
combined with an appropriate solver for the electric field, 
the algorithm is able to provide a full 'device solver' for 
a liquid crystal display. The method employs two lattice 
Boltzmann schemes, one for the evolution of the momen- 
tum and one for the evolution of the Q-tensor. This is 
necessary because of the large differences in time scale 
for the evolution of velocity and director fields in a typ- 
ical display or experimental arrangement. The paper is 
organised as follows. In Section^ we outline the Qian- 
Sheng equations for tensor nemato-dynamics. In Sec- 
tion ^H] we present an LB method to recover these equa- 
tions. In Section ITVI results are presented for the vafida- 
tion of the method against analytical equations and the 
application of the method to the modelling of a Zenithally 
Bistable Device (ZBD) device is reported. Section 
concludes by highlighting the benefits of the method de- 



scribed in the paper and discusses the implications for 
future work. A Chapman-Enskog analysis {eg 0) to 
justify the LB scheme is given in Appendix and Ap- 
pendix ^ summarises useful relationships between the 
vector, tensor and LB parameters of the method and lists 
the material constants used in the simulations detailed in 
Section CVI 



II. THE QIAN-SHENG FORMALISM 

In this section we summarise the Qian-Sheng formal- 
ism for the flow of a nematic liquid crystal with a 
variable scalar order parameter. The tensor summa- 
tion convention is assumed over repeated greek indices 
which represent three orthogonal Cartesian coordinates; 
no summation convention is assumed for roman indices 
which are used to indicate lattice directions of the LB 
algorithm. 5af3 and Sap-f are the Kronecker delta and 
Levi-Civita symbols respectively and a superposed dot 
( ) denotes the material time derivative: dt + Uada- 

The symmetric, and traceless, macroscopic order ten- 
sor, Q, is defined to be 
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where S and are the uniaxial and biaxial order pa- 
rameters with n, I and rh being orthogonal unit vectors 
associated with the principle axes of Q. In the uniaxial 
approximation P^ = and in the ELP approximation the 
scalar order parameter S Sq, a constant. The direc- 
tor, h = {sin 9 cos (j), sin 9 sin (j>, cos 9), is the eigenvector 
corresponding to the largest eigenvalue of Q. 

Following the momentum and order evolu- 

tion equations for incompressible {daUa = 0) nemato- 
dynamics are written as 



pup = da {-P6ap + <p + + aff ) 
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Here the local variables are p the liquid crystal density, 
u the fiuid velocity, P the pressure, J the moment of 
inertia. A and are Lagrange multipliers chosen to 
ensure that Q remains symmetric and traceless. cr^^ and 
hap are the distortion stress tensor and molecular field 
defined by the Landau-deGennes free energy, F, for the 
system through the expressions 
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a^p and h^^p are the viscous stress tensor and viscous 
molecular field respectively and are defined by 



(6) 



global free energy density is considered to be a sum of 
contributions arising from a number of different physical 
phenomena Tciobai = / Fsuik dr + J Fsurface dS. The 
free energy densities have the form 



Fsulk — FlcIG + FElastic + FElectr 
+ Fpi exo 



where 



(11) 



(7) 



Here (3i, fi^ are equivalent to the ELP viscosities, -/Va/3 is 
the co-rotational derivative, Nap — Qa/3 — s, 



eptiy^tiQay = \ {d^Up + dpUa) and Wap = 

\ {daUp — dpUa) are the symmetric and anti-symmetric 
velocity gradient tensors with the vorticity being lo^ = 
\s-yapWap- cr^p^ ^h*^ stress tensor arising from exter- 
nally applied electromagnetic fields 
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(8) 



where E (H) is the electric (magnetic) field strength, D 
the electric displacement vector and B the magnetic fiux 
density. 

Direct calculation of the trace and off-diagonal ele- 
ments of Eq. (O shows that lagrange multipliers are given 
by A = i {^hjj — ^^2^77) and = ^Sap-yhap- The term 
in is included in order to correct the small compress- 
ibility errors that arise in LB techniques when the condi- 
tion upon the Mach number (velocity to speed of sound 
ratio), M = \u\/cs ^ 1 is violated. 

Order of magnitude estimates and experiments both 
show the influence of the moment of inertia to be negligi- 
ble; we therefore set J = in Eq. Following J^, the 
viscous stress tensor and the equation of motion Eq. (O 
are recast in a form more suitable for the LB develop- 
ment. 
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The derivation of expressions for the molecular field 
and distortion stress tensor follows the phenomenological 
approaches for the free energy of liquid crystals ,27J . The 



FhdG — Fiso + -^apQapQpa — PpQal^Q (i-jQ-ja 

+lFQapQ0aQiiuQuii. (12) 
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p2 {.FiQ fivQvfi L^Q ^nQpT-Q-f (13) 

^ch 

'-^^Q^'^'^^ EaQapEp — -ea€j^E^ (14) 

FMagnetic = ^ -^fJ-OX'a"'^ FlaQapH p — -floXljFf^ (15) 

FfIcxo = —S,lEadjQaj — S,2EaQajdfiQj^ (16) 

Fsurface = -TT {Q a0 - Q°aaf (1^) 



The coefficients ap, f^FilF are parameters controlling 
the phase of the thermotropic liquid crystal, the negative 
(positive) sign preceding the [3f term dictates a calamatic 
(discotic) state; for biaxial phases sixth order terms are 
used. Li, i = 1...4, determine the elastic constants. 
Pch is the pitch of any chirality with /xg (eg) being the 
permeability(permittivity) of free space, x and e are the 
diamagnetic and dielectric tensors with x™°^ {^^"'^) the 
maximal diamagnetic (dielectric) anisotropy {ie S = 1). 
^1 and ^2 are flexoelectric constants, W an anchoring 
strength and a preferred surface state. This form for 
the free energy maintains consistency with the Q-tensor 
dynamics equations in that a direct analogy with the 
experimental ELP parameters is made (see Appendix IBl 
for the relation between experimental ELP values and 
the Q-tensor method). 

To close the governing equations at surfaces, non-slip 
boundary conditions are imposed upon the velocity. For 
infinitely strong anchoring a Q is specified according to 
Eq. . In cases of weak anchoring the order tensor at 
the surface evolves according to 



fJ-sdtQaP — h^p — X^ Sap — Eap-^X':^ 



(18) 



where h^p = 



OFbm j-, _ dFs^rfac. xS _ luS \S 
WTOlfi) ^ dQc0 ' ~ 3 "77' ^7 



^^a/i^h^p, is an outward pointing surface unit normal 
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vector and /is is the surface viscosity defined through 
A*S = A*i^S where Is is a characteristic surface length 
typically being in the range Is « 100 A 1000 A 28]. 

A non-dimensionalisation of the governing equations 
with respect to characteristic velocity, U, length , L, vis- 
cosity, rjeff — 5 — 3^^, and elastic constant, Zi, 

yields three key dimensionless numbers which govern the 
dynamics of the momentum, director, and order param- 
eter respectively 
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The characteristic timescales, f , for variations in the mo- 
mentum, director and order parameter are also given. Re 
and Er are the Reynolds and Ericksen numbers. De is 
the ratio of the relaxation time for the order parameter, 
pi/ap, to a time scale associated with the flow, L/U; it 
is similar to a Deborah number. Considering typical de- 
vice parameters, p ~ lO'^kgm"'^, rjeff ^ 10~^kg m-^ s~-^, 
L - IQ-^ni, tj ~ 10"^ms~\ Li - IQ-^^j^gms-^ and 
ap ~ lO^kgm^^s^^, we may estimate fp ~ 10^^ s, 
fn ^ 10~^ s, fs ^ 10~^ s. It is apparent the relaxation 
rate of the momentum compared to the director is much 
quicker, as is the relaxation of the order compared to the 
director and accounting for these timescale differences is 
essential for dynamic calculations. 



III. THE ALGORITHM 

We proceed now to describe the LB method which re- 
covers the set of equations set out in Sec. The algo- 
rithm is defined in Sec. IIII Al In Sec. IIII Bl we address 
the different timescales involved in liquid crystals' dy- 
namics and how to implement these in the LB method. 
A Chapman-Enskog multi-scale analysis of the LB algo- 
rithm is given in Appendix 1X1 



A. Statement Of The Algorithm 

LBGK algorithms {eg are well established for 

solving the Navier Stokes equations for isotropic flu- 
ids [301 ■ In order to recover the Qian-Sheng equations 
of Sec^lwe introduce two LBGK algorithms, one for the 
evolution of the momentum based on a scalar density 
fi {x,t) and a second LBGK scheme based on a tensor 
density giap (x, t) to recover the order tensor evolution. 
It is important distinguish between SI symbols in Sec. ^ 
and the symbols used in the LB algorithms, which are 
defined in terms of lattice units. However in this sec- 
tion, and Sec. Appendix IXI this distinction is ignored for 
clarity. In Sec. IIII Bl and Appendix IbI the distinction be- 
comes important and a prime is used to denote a lattice 



value. Further, a superscript ^ C^) is used to distinguish 
between momentum (order) algorithms. 

The principal reason for separating the momentum and 
order evolution algorithms is the very large difference in 
time scales between the two processes noted above. In 
each algorithm, forcing terms are used to recover the re- 
quired additional terms in the stress tensor and order 
evolution equations. This approach is more straight- 
forward to implement than the anisotropic scattering 
method used in an earlier work |2ll |. 

The LBGK algorithm for an isotropic fluid may be 
written in the form 

f,{x+c,At,t + At) = f,{x,t) (20) 

- ^ [n{x,t) - ft'H'^^t)) +^Ux,t) 

where fi {x, t) is the distribution function for particles 
with velocity Ci at position x and time and At is the 
time increment. Z/'^''' (a;, t) is the equilibrium distribution 
function and is the LBGK relaxation parameter. The 
algorithm fluid density and velocity are determined by 
the moments of the distribution function, 
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(21) 



The mesoscale equilibrium distribution function appro- 
priate to recover the correct hydrodynamics of incom- 
pressible fluids {M <^ 1) is. 



(eg) 
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where ti are lattice weights, ti, Ci, Cg are all dependant 
upon the choice of lattice, appropriate values of these 
parameters are summarised in |8(lj . An analysis of the 
standard isotropic algorithm identifies the lattice pres- 
sure and kinematic viscosity to be given by 



P 



2 



(2r^-l)At (23) 



(f>i is a forcing term which is chosen to recover the re- 
quired terms in the stress tensor. For a nematic liquid 
crystal governed by Eq. ||2Jl it is defined to be 



= UcixdpFpx 



(24) 



where 
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with analysis identifying (see Appendix lA l|l 



P = pc^ ' A^2A . 



i4 



, PC (2r,-l)At = /34-f^ (26) 



and a macroscopic observable velocity of pva = 
Yli fiCia + {^t/2)dpFf3a- The latter redefinition of the ve- 
locity is necessary to reduce higher order artifacts which 
are introduced by a position dependent forcing term '3l| ■ 
To recover the order evolution Eq. we retain the 
simple LBGK form but replace the scalar density fi {x, t) 
with a symmetric tensor distribution 3^/3 {x, t) evolving 
according to 

t) (27) 
— - (gtaiB {x, t) - gf^l (a;, t) \ + Xiafi [x, t) 

Here g^^^^ {x, t) is the equilibrium order distribution func- 
tion and the LBGK relaxation parameter for the order 
evolution. The lowest moment of the order distribution 
function, and its associated equilibrium function, are de- 
fined to recover the order tensor of unit trace. Sap 



'q/3 



E\ ^ (eq) 



(28) 



which is simply related to the dimensionless zero trace 
order parameter Q through the relation 



SSaf! — Saf3 



(29) 



The equilibrium order distribution is taken to be 



- UaUfj 



2c4 



(30) 



same lattice parameters defined for the 
momentum evolution. The forcing term Xiap is chosen to 
provide the rotational forces required correctly to recover 
Eq. ^ 



2ti At 
Xial3 = — 3— 
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The analysis (see Sec. lA 2\i identifies the key relation 



(31) 
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The scheme described here involves two coupled LB al- 
gorithms. Both may be run independently; for example 
if the effect of flow is to be ignored or only static equilib- 
rium configurations are desired, running the giap scheme 
alone will suffice. In practice for typical device geome- 
tries, the flow fields evolve on a much faster time scale 



than the director field; to model such systems the mo- 
mentum is evolved to steady state between each time 
step of the order evolution equation. Although the time 
taken for the momentum to reach equilibrium is signifi- 
cantly shorter than the time step of the order evolution 
equation, the loss of accuracy in this approach is small. 



B. Timescales In The Algorithm 



Constructing an analogous set of dimensionless num- 



from Eqs. 
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Re' 
Er' 

De' 


2U'pL'p 


- 1 2L'P^ 


- c2(2r^-l) 
2U'QL'Q 


- c2(2r^-l) 
^ ' - 2L'«" 




' " c^(2Tq-1) 
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(33) 







OS'' 


So 



We choose = = 1 and L'^ = L'« - 

The latter identity sets the simulation size to resolve 
variations in Q. The correct dynamics are achieved by 
matching the algorithm dimensionless numbers Eq. H33() 
to the real dimensionless numbers Eq. (|19|l . From 
Eq. 133|l this requires that C7'^ differs from U''^ by an 
amount Er / Re. In order to recover an internally consis- 
tent simulation, these different values of the LB velocities 
in the momentum and order evolution algorithms require 
the forces to be appropriately scaled when information is 
passed between the two algorithms. A list of the scaling 
is given in Appendix FBI 

We may now take the ratio of characteristic SI to LB 
times to give the time value of the LB discrete time 
step. Using typical values shows that: Ai^ ~ 10"^'^ 
s, IStn ^ 10^® s, Atg ~ 10~® s. Hence the momentum 
algorithm needs to be iterated many times within a single 
iteration of the order algorithm. Alternatively for lam- 
inar creeping fiows. Re <C 1, the equilibrium flow field 
will be reached in a small number of l\tp and we may 
jump forward in time to the next Ai„ ~ Atg reducing 
the overall processing time. 



IV. RESULTS 

In this section we begin by presenting results which val- 
idate the algorithm developed above by comparing its nu- 
merical predictions with analytical results for some sim- 
ple cases. We then show how the technique may be used 
to study the motion of defects in a commonly studied 
bistable liquid crystal device. 
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Comparison with analytical results for the 
Miesowicz viscosities 



We first consider ttie flow alignment of the director in 
a shear flow in the absence of an external aUgning field. 
Provided the channel width is sufficiently large, we may 
ignore gradients in Q in the centre of the channel. In this 
case the alignment at the centre of the channel is solely 
determined by the viscous torque. From Eqs.QJ and H1U|) 
we can solve for the director angle, 6, to find 



cos(26l) = -l^{iS + P^ 



71 f S+\P, 

72 V <S'o 



(34) 



A second standard case is to measure the shear viscos- 
ity of the nematic in the presence of a strong external field 
which imposes a fixed director angle. These experiments 
yield the Miesowicz viscosities 11]. However, the stan- 
dard results must be extended for the case of a variable 
order parameter. For an arbitrary fixed director angle 
the effective viscosity, 77*, is found to be 
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~ 1) - fS{inl - 1) 

ItLS^nlnl -f 2fi^2„2„2 _ 9|1 ^2^2„2 
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from which the following Miesowicz viscosities can be 
determined:- 



_ Pi 








2 


4 


4 




_ Pi 


1 3/i2S' 


P'oS 


I PeS 


2 


8 


4 


^ 2 


_ /34 






PeS 


2 


8 


4 


2 



at fia = (0, 1, 0) 




Vb 

Vc 



these being identical to the EL expressions in the 
limit S* — > 5o. A biaxial correction is not required as the 
aligning field serves to cancel biaxial contributions from 
the shear. 

In order to assess the accuracy of the method described 
in Section IIIII we tested it against these analytical val- 
ues. We used a channel width L = 1.2/im, a shear rate 
7 = 10'*s~^, viscosities {ai — —0.011, a2 — —0.102, — 
-0.005, a4 = 0.074, as = 0.084, ag = -0.023} kg 
m^^ s^^. Landau parameters { a = 65000 J m~^ K~^, 
B = 530000 J m-3, C = 980000 J m'^} and T = 
Tin — 4(T/jv — T*). The boundaries were assumed to 
have infinite anchoring and the flow induced by adding 
2Up'"c^-2£^/{cl) to the right hand side of Eq. ^ where 
the wall velocity is +{—)u^ at the top (bottom) bound- 
aries with periodicity in the x and y directions. 

In the absence of an aligning field, the director angle 
in the shear flow was found to be 12.166° which agreed 
with the value predicted by Eq. to 7 signiflcant 

figures. Accuracy was found to be maintained over all 
flow aligning viscosity ratio's with a typical increase in 
S around 0.002 and biaxiality = 0.002. The Miesow- 
icz viscosities were measured using an aligning field of 75 
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FIG. 1: Schematic of the two dimensional ZBD geometry over 
two grating pitches, w. Homeotropic boundary conditions 
serve to cause bistability. Simulations contain one grating 
pitch and periodic boundaries in the x and y directions. 



volts (Aea = 10.3) in the relevant directions. Non-slip 
boundary conditions were applied using the bounce-back 
method [s^l) and the fiow was induced by applying a 
constant body force at z = L/2. The resultant viscosity 
ratio's are compared in Table^where data is measured at 
z = L/A. It can be seen that the LB solver gives results 
in good agreement with the expected values. 



Theory 


Simulation 


% error 


riahb = 1-446 


7a-V76"' = 1-446 


1.6 X 10"" 


rial lie = 0.227 


7a"V7c-' = 0.227 


1.4 X 10"" 


r^b/ryc = 0.157 


75'V7c"' = 0.157 


2.8 X lO"'^ 



TABLE I: Table of theoretical (Eq. I|36|l 'l and simulated ratio's 
of the Miesowicz viscosities. 



Investigation of defect motion in a bistable 
device 



B. 



The ZBD device [2|32| uses a structured surface, such 
as that in Fig. to introduce bistability which may 
be used to design a very low power display. The two 
bistable states are characterised by the presence or ab- 
sence of defects which will be referred to as the defect 
(2?) and continuous (C) states respectively. One possi- 
ble method to latch switch between states is to use an 
electric field that couples to the fiexoelectric properties 
of the liquid crystal material. Latching between the two 
states is a dynamic process which involves the nucleation 
and annihilation of defects. 

We followed |2| and modelled the ZBD surface with 
the function g{x) = /i sin -|- ^ sin (2^)] projected 
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FIG. 2: Defect trajectories during the V to C latching for various En values. Points indicate location of the annihilation. 
The inset figure indicates the time at which defects annihilate from the turn on of the voltage, e^-, — 18, V = +18 volts and 
Aca = 10.3. 



onto the LB boundary over one grating period, w; the 
height of the grating is h and the parameter A controls 
the level of asymmetry. Weak anchoring conditions was 
used at the surfaces by implementing Eq. H18|l with an 
explicit forward time finite difference method. The gra- 
dients in the equation were extrapolated to second order 
from the bulk and an average taken over the values ob- 
tained from each lattice direction. It should be noted 
that in order to achieve equality of the elastic constants 
in the bulk and at the surface it is necessary for the sur- 
face parameters Lf to be different from bulk LB through 
theLf = Li/(2rQ-l). This arises because the relaxation 
processes in the bulk, which are governed by the param- 
eter Tq , contribute to the measured elastic constants in 
the bulk. However, there is no equivalent collision pro- 
cess in the surface algorithm. 

In the presence of the voltage applied to the device, 
it is necessary to solve Maxwell's equations over the LB 
grid to obtain the local values for the electric field, E. 
For completeness these equations are 

dcD^ = 1 

Da = f'Q^apEp + Pa I /o7\ 

£a0 = (SAe^^'^Qa/J -t- e.yjSal3)/3 _ 

in which V is the local voltage, e-y^ — 2e± + e||, Ae^ = 
S'Ae™"^ and is defined from Eq. (|16|l by writing it in 
the form Fpiexo = —PaEa. We solve equations (|S7|l using 
a successive overrelaxation method at each iteration of 
the LB algorithm for giap- This therefore determines the 
electric field which is consistent with the instantaneous 
value of the Q tensor. 



We investigated the effect of material properties on the 
motion of defects in this device; in particular we studied 
the interplay of dielectric, flexoelectric and surface polar- 
isation effects. We used the set of material parameters 
given in Appendix ^ The system was first established 
at steady state in one of the equilibrium states; the simu- 
lation was then run using the algorithm described above. 
The equilibrium states were located by starting from an 
appropriate initial condition and running only the giafi 
algorithm. The defect equilibrium state {T>) has a -1/2 
defect near the peak of the grating and a -1-1/2 defect 
near the trough of the grating. 

In Fig. © the flexoelectric coefficient S13 = £ii±£m , 
(eii — 633) is varied. Starting in the V state and applying 
-f 18(0) volts to the upper (lower) electrodes the resultant 
defect trajectories are shown in the grating region. For 
i?i3 = the defects move slowly along the surface and 
annihilate. As i?i3 increases we increase the surface po- 
larisation and order parameter which pushes the defects 
further out into the bulk of the device to annihilate. As 
i?i3 increases the — i defect mobility is increased as seen 
in the annihilation locations. The inset of Fig. Q shows 
the time taken for the defects to annihilate which is an 
indicator of the latching speed. It is found increasing Ei^ 
increases the latching speed. 

Alternatively we may keep -E13 constant and vary the 
dielectric anisotropy, see Fig. For a deceasing Aeq we 
effectively increase flexoelectric contributions to the ne- 
matic (see Eqs. l(Tl|) and this increases the surface 

polarisation that pushes the defects further away from 
the surface. Increasing Ae^, effectively reduces the flexo- 
electric contributions and the defects annihilate closer to 
the surface; this also reduces the latching time. 
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FIG. 3: Defect trajectories during the D to C latcliing for various Aea values. Points indicate location of the annihilation. 
The inset figure indicates the time at which defects annihilate from the turn on of the voltage. e~f^, = 18, = +18 volts and 
Ei3 = 20pC m-\ 
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FIG. 4: Defect trajectories during the T> to C latching for various e^^ values. Points indicate location of the annihilation. The 
inset figure indicates the time at which defects annihilate from the turn on of the voltage. Aea = 10.3, V = +18 volts and 
£13 = 20pC m-\ 



Fig. (01 has fixed Ae^ and Ei^ but the grating permit- 
tivity e;^^ is changed. This has the effect of diffracting 
the electric field lines for an increased mis-match of sur- 
face and nematic permittivities. At the lower dielectric 
mismatch the defect annihilation location is at ~ h/2. 
Increasing the dielectric mismatch increases the mobility 
of the — i defect allowing it to travel further and anni- 
hilate near the grating trough. There appears to be an 



optimum value of e;^^ ~ 26 for which the annihilation 
time is shortest. 

Fig. ||SJ) shows the effect of increasing the applied volt- 
age for constant e^^, Aea and E13 ; this has the effect 
of increasing the contributions of both the flexoelectric 
and dielectric terms. An increased voltage tends to cause 
the defect trajectories to move away from the surface to- 
wards a saturation distance for which further increase 
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FIG. 5: Defect trajectories during the to C latching for various E13 values. Points indicate location of the annihilation. 
The inset figure indicates the time at which defects annihilate from the turn on of the voltage. Aea = 10.3, e^~f — 18 and 
Ei3 = 20pC ia-\ 



causes little difference. Above the voltages shown in the Acknowledgments 
figure, a different latching mode is seen in which several 
pairs of defects occur in the annihilation process. As with 
a Fredericksz response an increased voltage results in a 

faster latching response. We thank Dr. I.HaUiday, Dr. D.J. Cleaver, Dr. S.V. 

Lishchuk, from Sheffield Hallam University and Dr. J.C. 
Jones, Dr. R. Amos from ZBD Displays Ltd for very 
useful discussions and comments in the developments and 
investigations undertaken in this paper. 



V. CONCLUSIONS 



An LB method has been presented which can be used 
to predict the dynamics of a nematic liquid crystal in the 
complex geometries which are increasingly being adopted 
for display devices. Nematic order, director, velocity, 
electric fields and surface polarisations are all recov- 
ered; this allows comparison to experimental results to 
be made for a wide range of cell geometries or surface 
patterning, fn the presence of structured surfaces and 
defects, it is essential to consider the variation in the or- 
der parameter. Essentially, a full 'device solver' has been 
developed and example results are given that show both 
the accuracy of the solver and its use in determining be- 
haviour of next generation LC devices. The influence of 
surface polarisations resulting from dielectric and flexo- 
electric effects are shown to effect defect trajectories and 
ultimately latching speeds. The solver is currently be- 
ing used in the development of next generation bistable 
devices. 



APPENDIX A: CHAPMAN-ENSKOG ANALYSIS 
OF THE LB ALGORITHM 



In this section we present a Chapman-Enskog analysis 
of the momentum and order evolution schemes. This 
analysis serves two purposes: to demonstrate that the 
method recovers the required governing equations and 
to identify the relation of the LBGK parameters to the 
associated transport coefficients and forcing terms. 
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1. Momentum Evolution 

The moments of the distribution function, fi, are de- 
fined to be 
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The the velocity basis and the are chosen to give 
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(A2) 



where 



(A3) 



Using a Taylor expansion on the left hand side of Eq. 120|) 
we obtain: 



(A4) 

We assume a forcing term of the form 0^ = tiCi\dpFp\ 
and use a multi-scale expansion, to second order 



h ^ st , t2 ^ e^t , dt = edti + e^dt^ 



xi — ex 



(A5) 



Using this expansion in Eq. I|A4|I and collecting terms we 

obtain: 

O(e0) 



j(0) ^ j(eq) 



(A6) 



(A7) 



1 

T 

2 



(1) 



.At 



At {dt, + Ciadajf, 
-Tp At [dt^ + Ciada^ ) °^ 

+TptiCixdp^Fpx 



p(2) 



(A8) 



in which we have used the 0{e^) result of Eq. IjATp to re- 
place a term of the form {dt, +Ciada, in the 0(£^) re- 
sult. Taking the zeroth moment of Eq. I|A7|I and Eq. (|A8|I 
whilst respecting Eq. lAip yields: 



0{e^) dt,p + da, {pua) ^ 



(A9) 



0{e'^) dt^p + 9^2 {pua) = cldj,d0,Fp^/2 
which can be recombined to give the continuity equation: 



dtp + da {pUa) = 



(AlO) 



where the term cld\dpFf3\/2 is corrected for by redefin- 
ing the macroscopic velocity as described below Eq. H2t)|l . 

Taking the first moment Cifj) of the first and 
second order Eq. (|A7|) and Eq. ljA8|) whilst respecting 
Eq. I|A1|I yields: 

0{e^) dt, {pup) + da,Tl^al = ^^71^7/3 

0{e^) [l-^)da,lV% + dtApup) 

-\-da2^\Jp = ■^4^72^7/3 fdt,dj,Fjp 

(All) 

In order to progress, the Tl[^p term needs evaluating and 

this requires knowledge of f^^'' (in Eq. (|A7|) ) . Using 
Eq. to 0{u) in Eq. (|A7I) . taking its zeroth moment 
and back substituting the result {dt,p = —dp,{pup)), fol- 
lowed by taking the first moment and another back sub- 
stitution {dt,{pup) = —c1dp,p — Tpc1d~f,Fjp) yields 

„(1) _ TpAtdp,{pUa)Hiap 
Ji ^ 2 

+ Tpt,c,x{TpAt + l)d^,F^x (A12) 
in which the symmetric quantity Hiap is given by 



FliaP = tiC 



.2 



SaP 



(A13) 



The symmetry of Eq. ljA13p allowing us to replace 
d\{pup) by pA\p in Eq. (jA12|l . in the incompressible 
limit. 

We now use Eq. (|A12|I to evaluate Eqs. HA11() and com- 
bine the results to find 

dt{pup) + Uada{pup) -dp{pcl) 

+c2(2Tp - l)Atda{pApa) + ^^^aFap (A14) 

Here the term —-^dtd-yFpj can be neglected assuming 
high order gradients are negligibly small. 

A detailed comparison of the terms in this equation to 
the target momentum equation Eq. Q gives the identifi- 
cations made in Eqs. H2()l) . Note that the isotropic terms 
may be incorporated into the scalar pressure |l ij . Com- 
parison of the remaining stress tensor terms allows the 
forcing term Eq. (|25|l to be defined. This completes the 
momentum analysis. 
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Order Tensor Evolution 



which may be written in terms of Q as 



The moments of the order distribution function are 
defined to be 



1 




Sa[3 
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_ "q/37 _ 



n > 



(A15) 



It should be noted that we adopt a unit trace order tensor 
in these moment definitions in order to be consistent with 
|2ll |. However, since the momentum and order are now in 
separate algorithms, we could equally well have defined 
zero trace moment definitions as in 20]. 

Using a Taylor expansion on the left hand side of lattice 
evolution equation (Eq.((23)) we obtain: 



At 



Q V 



2 



(A16) 



We suppose the as yet unknown forcing term Xi/^i/ will 
be dependent on the gradient in both u and Q and can 
be expanded as Xitiv 
Eqs. (|A5|) with the expansion 



^XiX + s^xfX- We augment 



9^ 



:„W +£2 (2) 



(A17) 



Substituting into the Taylor expansion Eq. (|A16|I . we find 



(0) ^ (eg) 



(A18) 



9^ (A19) 



(1) (2) _ (2) .. . 

in which we have used the O(e^) result of Eq. I|A19|I to 
replace a term of the form [dt^ +(^iadai)9i^,y in the O(e^) 
result. 

Taking the zcroth moment of the first order Eq. (jA19|) 
expansion and using Eq. l|IT5)l gives 



(1) 



^ At 



(A21) 



Similarly, the zeroth moment of the second order expan- 
sion gives 



1 \ 2 2 

+ iEx£l-^Ex!;.l (A23) 



We now need to evaluate ^^al,v by obtaining an expres- 
sion for in Eq. (|A19|) . We use Eq. to 0{u) in 
Ea. (|A19|l . Taking the first moment of this gives 

i 

^u^dt,{S^,)^d^,{S^,cl)\ (A24) 

Using the result Eq. HA2ip . we may replace dt,{S^^) and 
find 



-T,^At 



+^qEc'7x!;.1 (A25) 



We may use the result obtained in text above Eq. (jA12|) 
to find: 



OtiUp — (AZb) 



and hence from Eq. ljA25|) and the incompressibility con- 
dition we find 



E + ^'d,. is,.)] + r,, ^^7X£1 (A27) 

i / i 



Upon converting 5 to Q Eq. I|A27|I is inserted in the 
earlier second order zeroth moment Eq. ljA23p giving: 
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TpT AtSf,^cl rp s , 2r / a \ a / (l)^ 

^ Oai[Op^fl3a) H Oai(UaU^ 

2T„Atc2 



(A28) 



\ ^ (1) \ ^ 



In order to simplify Eq. IA28I we omit terms which 
include third order gradients in either u or Q. Further 
in the limit M — ^ I, which holds for low Re LCs, 
we may omit the term which includes the product u u. 

Recombining the 0{e^) Eq. (|X22|l and ©(e^) Eq. (|X28|) 
expansion we obtain: 

c2 

dtQ^^ + UadaQ^,v = y (2tq - l) Mda (<9qQ^,.) 

~ Xiuv ~^ 2/\t '^""^ (A29) 

A comparison of the terms in this equation and the 
target order equation (Eq. gives the identification 

made in Eq. H32|l . We now compare the remaining terms 



with the force terms in Eq. (|A29p . We make the as- 
sumption that the forcing term must be introduced at 
O(e^) since it is gradient dependent; we therefore choose 
X^i Xi^u = 0- This allows us to make the identification 
for the forcing term, Xiap^ given in Eq. 

APPENDIX B: ALGORITHM PARAMETERS 

Here we give details on the relations between EL ma- 
terial coefficients and the Q tensor coefficients. We use 
standard EL notations as given from [l3|. The details 
are obtained by using Eq. ^ in Eqs. © and ©. Note 
So stands for the equilibrium order parameter, not sim- 
ulation evolved order parameter, S. 



_ 2(03-02) _ 271 2(02-1-0:3) _ 272 — a _ n II — 

/^l ~ 9Sl ~ 9Sj ' — 3Sj, — — l~>Q 1-^5 , f^s — 957 

Ll = 2^ {3K22 + — Kii) , L2 = gl^- {Kii — K22 - K24) , L3 = ■^K24 , L4 = 2^ (^33 — Kii) 

ap = — T*) , Pp — , -fp = |C , Sjn — ^ , T* — Tjn — 

^1 = 91^(^11+2633), 6=91^(611-633) 

(Bl) 

Then the relation of the Q tensor coefficients to both momentum algorithm and the order algorithms are 
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Note that in the definition of Atg we have used an elas- 
tic constant which is characteristic of a simple twisted 
nematic cell to illustrate the mapping of variables, plus 
eg" = 1 and A/^ = At''' = At' and L'" = L'" = L' and 

^P^i^^JMf and</>^0^Qj ^^^^^ 

The material parameters used for simulations in tJlV B\ 
are: {K^ = 10, ^^22 - 7, ifgs = 14, 1^24 = 5) x 
IQ-i^kg m s"^, (ai = -11, = -102, ag = -5, q;4 
74, as = 84, = -23) x lO^^kg -1 s-^, p = 



1.01 X lO^kg m"^ T ^ Tin - 4:{Tin - T*) K, a = 
65000Jm-3K-i, B = 530000 J m^^^ C = 980000 J m-^ 
, W = 7xl0-^kgs-^, ea = 10.3, e^^ = 18, 5(611+633) = 
2x 10^^^^ S" m~^, ^5 = lG^''m. Other specific constants 
are provided in the figure captions. The parameters used 
represent a hybrid of commonly used materials; they do 
not correspond to a specific material since a complete set 
of material parameters does not exist in the literature for 
one material. 
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